Nonequilibrium Green's function method for quantum thermal transport 
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This review deals with the nonequilibrium Green's function (NEGF) method applied to the prob- 
lems of energy transport due to atomic vibrations (phonons) , primarily for small junction systems. 
We present a pedagogical introduction to the subject, deriving some of the well-known results such 
as the Laudauer-like formula for heat current in ballistic systems. The main aim of the review is 
to build the machinery of the method so that it can be applied to other situations, which are not 
directly treated here. In addition to the above, we consider a number of applications of NEGF, not 
in routine model system calculations, but in a few new aspects showing the power and usefulness 
of the formalism. In particular, we discuss the problems of multiple leads, coupled left-right-lead 
system, and system without a center. We also apply the method to the problem of full counting 
statistics. In the case of nonlinear systems, we make general comments on the thermal expansion 
effect, phonon relaxation time, and a certain class of mean-field approximations. Lastly, we exam- 
ine the relationship between NEGF, reduced density matrix, and master equation approaches to 
thermal transport. 
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I. INTRODUCTION 

The method of nonequilibrium Green's functions 
(NEGF) was initiated by Schwinger in a rather math- 
ematical paper [l| for a treatment of Brownian motion 
of a quantum oscillator. Already in 1961, the impor- 
tance of forward and backward evolution in the calcu- 
lation of nonequilibrium quantum expectation values at 
time t evolved from an earlier time was recognized and 
the six different Green's functions defined. The next im- 
portant development in NEGF came due to Kadanoff and 
Baym Q where the main emphasis was to derive quan- 
tum kinetic equations. Keldysh showed that diagram- 
matic expansion is possible even for nonequilibrium pro- 
cesses [3|], a key idea being contour order. These initial 
developments all occurred in the early 1960s. There are 
a number of earlier reviews 0-111 and conference series 
0, [1] that people working in this field should be aware 
of. An important paper on treating transport by NEGF 
is that of Caroli, et al. where for the first time, an 
explicit formula for the transmission coefficient in terms 
of the Green's functions was given. Its modern form pre- 
sented here is due to Meir and Wingreen [l(J|. Some of 
the very recent reviews on NEGF method, mostly still 
for electronic transport, can be found in Refs. [ll| - [l4l . 

This paper can be thought as an update to our earlier 
review [15] on the application of NEGF to phonon trans- 
port. The main aim is to develop the theory more sys- 
tematically and to review the various new applications. 
Some of the straightforward, routine recent applications, 
e.g., Refs. [I6I - I21I will not be discussed. We start slow 
with the problem of the harmonic oscillator in Sec. IL 
Since any phononic systems at the ballistic level can be 
thought of as coupled oscillators, and in eigenmodes, in- 



dependent oscillators, the single mode oscillator is funda- 
mental to the NEGF method. In Sec. Ill, we define the 
nonequilibrium Green's function proper. Here, we look at 
the contour ordered Green's function as well as operators 
used to define it by carefully introducing a contour ver- 
sion of the evolution operator as well as giving a formal 
definition of the Heisenberg operator on the contour. The 
mathematical aspect of functions defined on the contour 
is dealt with in Sec. IV. Two methods are available for 
obtaining equations and actual computation, the equa- 
tion of motion method and Feynman diagrammatic ex- 
pansion. Both of them are formulated on the contour of 
a finite segment [toj^M]- This is discussed in the follow- 
ing two sections, V and VI. We emphasize the view that 
the contours are defined on finite segment. This point of 
view makes the theory valid both for transient and steady 
state. The current formulas are derived in Sec. VII. The 
remaining sections review some applications, including 
problems of multiple leads, full counting statistics, which 
is to look at the full distributions of transferred energy 
in a given time interval. We review few applications in 
nonlinear situations where NEGF gives reasonably good 
results, this includes thermal expansion and phonon life 
time, and a self-consistent mean-field theory for a quartic 
nonlinear junction. NEGF is normally concerned with 
Green's functions, but it can also say much on the re- 
duced density matrix; here in Sec. XII, we review Dhar, 
Saito, and Hanggi's method of computing the reduced 
density matrix in steady state for a transport system. 
This is quite relevant with respect to the last topic of 
this review, the quantum master equation approach. We 
try to rephrase the usual quantum master equation in 
terms of NEGF and offer an approach and formula to 
obtain higher order current with respect to the system- 
bath couplings. We end the review with a brief summary 
in Sec. XIV. 
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II. HARMONIC OSCILLATOR, EQUILIBRIUM 
GREEN'S FUNCTIONS 



In this review, we take a bottom up approach to 'build' 
the nonequihbrium Green's functions from the equihb- 
rium ones. This will be done in the first few sections. 
In this section, we review the basic properties of a single 
degree of freedom harmonic oscillator in thermal equi- 
librium. The set of functions defined here will be found 
of great utility later as any phononic system (and even 
photonic systems) can be thought of as a collection of 
independent harmonic oscillators if we work in the eigen- 
modes. Hence the problem of the harmonic oscillator is 
fundamental to phonon transport. Similar discussions 
and formulae can also be found in the lecture notes 
of Brouwer Sec. 2.6 (an excellent introduction to 
NEGF), in Kleinert Chap. 18.5.1 (although the con- 
ventions are different from ours), as well as in appendix B 
of Ref. ^24: and in appendix A of Ref. 25. 

We assume that the reader is familiar with the solu- 
tion of the quantum-mechanical problem of a harmonic 
oscillator using creation/annihilation operators (see, e.g., 
Bohm |2^, Chap. II. 3). The Hamiltonian of a single 
quantum oscillator is given by 



H = 



P 

2m 



1. 



(1) 



where x is the displacement operator, p is the conju- 
gate momentum such that [x,p\ = xp — px = ifi, h is 
the reduced Planck constant, m is mass, and k is the 
force constant. For notational simplicity, in this review, 
we'll always perform the transform u — X\fm, so that 
the mass can be transformed away. Also introducing 
k = toJI^, where f2 > is the oscillator angular fre- 
quency, the Hamiltonian can be rewritten as 



1 



1 



1 



(2) 



where the mass-normalized displacement u can be ex- 
pressed in terms of the annihilation and creation opera- 
tors as 



2^) + 



(3) 



We have the commutation relation [a, a^] = 1. In the 
Heisenberg picture, the operators evolve in time, and the 
states do not change. The Heisenberg equation of motion 
for a takes a very simple form 



da{t) 
dt 



1 



'a{t),H] = -ina{t). 



(4) 



This gives the oscillatory solution a{t) = ae~*^*, thus 
the Heisenberg solution for u can be easily obtained. 

In equilibrium statistical mechanics, we assume that 
the system is not in a pure quantum state, but in various 
states with some probabilities. More precisely, we need 



to describe the system with a density matrix (see Huang 
[27} . Chap. 8). In this review, we'll always use the canoni- 
cal density operator, p = e'^" /Tiie-'^"), /3 = l/ikeT), 
where kB is the Boltzmann constant and T is the ab- 
solute temperature. In the energy eigen-state represen- 
tation, \n), the Hamiltonian is diagonal, and using the 
facts 



H\n) = {n+]^nVl\n), 

a\n) — ^/n\n — l), 
a^\n) = y/n + l|?i + l), 



(5) 
(6) 
(7) 

(8) 
(9) 

where the angular brackets denote trace with the canon- 
ical density matrix, i.e., (• • •) = Tr(p- • •), and 



we find 



(aa) = 0, (a^'a^) = 0, 
(at a) = /, (aat) = 1 + /, 



/ 



1 



(10) 



is the Bose-Einstein distribution function. 

Now we are ready to define correlation func- 
tions, or Green's functions for the harmonic oscillator. 
One can define the Green's functions using the cre- 
ation/annihilation operators — this is traditionally done 
in many-body theory. But for phononic systems, it is 
more efhcient if we only use the displacement operators. 
We define the greater Green's function as 



g>{t,t')^~-{uit)uit')), 



(11) 



where the time-dependence is from the Heisenberg evolu- 
tion, and the angular bracket is for the average over the 
equilibrium density matrix p. Using the relation between 
u and a, Eq. ([3]), and the solution of Eq. (|4]), and Eq. ([8]), 
(El), we get 

^^^^'^'^ " [/e^"(*-0 + (1 + /)e-^"(*-*')] . (12) 

We note that the function g^{t,t') — g^{t — t') is ac- 
tually a function of one argument due to time transla- 
tional invariance. It is always so provided that the system 
is in thermal equilibrium or in a nonequihbrium steady 
state. The greater Green's function is nothing but the 
position-position correlation function in time. The extra 
factor (—i/h) is there to match the evolution operator 
in a Dyson expansion and is purely conventional. With 
this definition, g^{t) has the dimension of time. We are 
going to define a few more functions: the lesser Green's 
function 

g<{t,t')^~^{u{t')u{t))=g>{t',t), (13) 
the time-ordered (causal) Green's function 
g'{tX} = -l{Tuit)u{t')) 

= 0{t - t')g>{t,t') + 0{t' - t)g<{t,t'), (14) 
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and the anti-time-ordered Green's function 

= e{t' - t)g> {t, t') + e{t - t')g< {t, t'). (15) 

The Heaviside step function is defined a.s 9{t) — 1 ii t > 
and otherwise. The meaning of the time-order operator 
T and anti-time-order operator T is given by the second 
fine. For the time order, the positions of the two oper- 
ators are unspecified until the time t and t' are known. 
The position of the operators is such that the operator 
on the right is the earliest and following the order of time 
as one goes from right to left. The anti-time order is the 
opposite. 

It seems redundant at this point to introduce four of 
these functions, as Eq. ([T^ determines all the others. 
However, these four functions form the components of so- 
called contour ordered Green's function, (;(r, r'), which 
has great utility when we deal with nonequilibrium situ- 
ations. Another pair of important Green's functions are 
the retarded Green's function, given by 



the domain from —oo to -I-cxd. The Fourier transform of 
the retarded Green's function for a single oscillator is 



9'M 



+ 00 



1 



{lo + iviY 



{t] -> o+). 



(21) 



It is important to add a small positive damping factor 
T] so that the integral converges. This choice displaces 
the poles in the complex plane of frequency below the 
real axis and produces the desired causality property that 
g{t) = for t < when one performs an inverse Fourier 
transform using contour integral. The advanced version 
is obtained by complex conjugate, (/"[w] = ^''[w]*. 
The lesser Green's function in frequency domain is 

7 TT 

5< H = - If^i^ + + moj + n)] . (22) 



Using the Plemelj formula 



1 



P- 



iiTS{x), 



(23) 



n 



and the advanced Green's function 



(16) 



g^{t,t') = -e{t'-t){[u{t),u{t')]). (17) 

The retarded Green's function appears in linear response 
theory, and it has the same meaning as that of Green's 
function in classical physics, i.e., it is the solution of the 
equation 



where P stands for Cauchy principal value, we can 
now relate the lesser Green's function with the retarded 
Green's function as 



3<M 



(9' 



9 



(24) 



This equation turns out to be true in general (in the 
sense of the G defined in the next section) for equilib- 
rium systems and is one particular form of a fiuctuation- 
dissipation relation. Another, which is actually equiv- 
alent to that one, is g^[a;] = e^^'^ g^[uj\. In the time 
domain, this is the Kubo-Martin-Schwinger (KMS) con- 
dition [gllli 5<(i) = g<{~t + ipn). 



■f{t) + n'g^t)^-S{t), 



(18) 



with the condition g^'{t) = for i < 0, where 6{t) is the 
Dirac S function. 

In practical calculation, particularly in the case of 
time-translationally invariant situation, it is more con- 
venient to work in the frequency domain. We thus define 
the Fourier transform of the Green's functions, e.g., by 



9'M 



and its inverse 



+ CXD 



(19) 



(20) 



We use the same symbol for a function of time defined in 
the whole domain (— cx), -t-oo) and its Fourier transform. 
Whether it is the function of time or its Fourier transform 
is indicated by its argument (t) or [to] . Just like the time 
t, since w is a Fourier transform variable, it also varies in 



III. NONEQUILIBRIUM GREEN'S 
FUNCTIONS, BASIC DEFINITIONS AND 
PROPERTIES 

In this section, we generalize the definitions for the sin- 
gle degree harmonic oscillator and consider a general sys- 
tem described by vibrational displacement Uj , where the 
single index j runs over all the relevant degrees of free- 
dom of the problem. For example, in a three-dimensional 
system, j may refer to the ^-th atom in the x direction. 
This compact notation makes the formula valid for any 
dimensions. We define the greater Green's function 
as a matrix, with the elements 



(25) 



where the trace is the quantum-mechanical trace over a 
complete set of states, Uj{t) is the Heisenberg operator 



for the displacement given by 

Uj{t) 



(26) 
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where Uj is the associated Schrodinger operator, and H 
is the Hamiltonian of the system. If the Haniiltonian is 
exphcitly time-dependent, one needs to replace the expo- 
nential factor by the full Schrodinger evolution operator 



U{t,to) = Texp 



Hf,dt' 



{t > to), (27) 



i.e., Uj{t) = U{tQ,t)ujU(t,to). Anti-time order needs to 
be used if t <tn in the above formula. We refer to Fetter 
and Walecka [30|, Chap. 3, for an excellent exposition for 
the three pictures in quantum mechanics and the prop- 
erties of the evolution operators. The density matrix 
in Eq. (I25p is at time to. Since p{to) here is arbitrary, 
the system in general is not in equilibrium, and the two- 
argument function depends on the two times, t and t' , 
separately. By 'nonequilibrium', we'll simply mean that 
p{to) is not given by a canonical distribution, cx e~^^ , or 
the Hamiltonian defining the dynamics may be explicitly 
time-dependent. Note that a reference time, to, when the 
Heisenberg picture and Schrodinger picture synchronizes, 
is arbitrary. Common choices are either setting to to 0, 
or the limit to — oo. 

Other Green's functions are defined in a similar fash- 
ion. The lesser Green's function can be obtained by swap- 
ping time arguments and space indices simultaneously, 

G<fc(i,0 = G>(t',t), (28) 

and the retarded Green's function is obtained by the com- 
mutator, 

q,(t,i') = -^0(t-t')Tr{p{to)[u,{t),u,{t')]} 



9{t-t') G>{t,t')~G<(t,t') 



(29) 



Similarly, advanced Green's function can be obtained by 
swapping arguments: 



G°jk{t, t') — Glj{t' , t). 



(30) 



The time-ordered and anti-time-ordered Green's func- 
tions can be obtained from the others defined above, 
= G< + and G* = G< - G°, as matrix equa- 
tions. We collect some of the useful relations among the 
Green's functions. 



G' + G^ = G> +G<, 
G'' ~ G" G> - G<, 
G'' + G° = G* - G*. 



(31) 
(32) 
(33) 



These linear relations are valid in both the time domain 
and frequency domain. In addition, an important rela- 
tion in the frequency domain is 



G° M ^G'^M^ 



(34) 



where the dagger f stands for hermitian conjugate. The 
retarded Green's function is analytic in the upper half 
plane of the complex frequency domain. This property 

















T 
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FIG. 1. Contour C used to define the nonequifibrium Green's 
functions. The upper branch is called -I- and lower one — . 
The order follows the direction of the arrows. 



guarantees the Kramers-Kronig relation relating the real 
part of G'"[cl'] with the imaginary part of G'"[ti;] (or vice 
versa) through Cauchy principle value integrals (we refer 
to Kubo et al. [31], Chap. 3.6, or Atland and Simons 
[131, Chap. 7). 

We have not yet finished with our definitions of Green's 
functions. The last and perhaps the most important 
Green's function in NEGF is the contour-ordered Green's 
function. The contour-ordered Green's functions are ex- 
plained in some books on many-body physics, e.g., Haug 
and Jauho [33|, Chap. 4, Zagoskin [3^], Chap. 3.4, Ram- 
mer [Hi, Kleinert [H, Chap. 18, or Di Ventra [11], 
Chap. 4, and Kamenev js^]. The usefulness of this type 
of Green's functions is because quantum evolution (the 
Heisenberg operators and density matrices) is two-sided, 
see Eq. , where we can think of Uj {t) as developing 
from the reference time to to the time of interest, t, using 
U{t,to), meeting the Schrodinger operator Uj at time i, 
and then being evolved backward in time by U{to, t) from 
t to to- The evolution goes forward and backward, form- 
ing a loop, or contour. Another deep reason is that only 
in this form of contour order, we can develop a transpar- 
ent perturbation theory, using the interaction picture. 

By convention, we define the contour G as going from 
to in the upper branch (forward going, -I-) or slightly 
above the real axis on a complex time plane, up to a 
maximum time tA[ relevant to the problem, then return- 
ing to the original time to from the lower branch (back- 
ward evolving, — ), or slightly below the real axis, see 
Fig. [TJ However, the time t is always real, and this has 
nothing to do with analytic continuation. We'll use the 
Greek letter r to denote a particular point on the con- 
tour and it is equivalent to a time t and a branch index 
a — ±. A nice, compact notation for this is f . An evo- 
lution operator U can be defined on the contour. If both 
Ti and T2 are on the upper branch with ri precedes T2 
(i.e. ij^ < ij), we define C/(t2, Ti) = C/(i2, ^i)- This is the 
usual evolution from a time ti to a later time t2 with time 
ordering [Eq. (P7)) ]. If both ti and T2 are on the lower 
branch with ri precedes T2, then t^ > ^2^, and U{t2,ti) 
is given by the same formula U{t2,ti), but since t2 < ti 
the time order should be replaced by anti-time order. If 
Ti is on the upper branch and T2 on the lower branch, 
we define U{t2,ti) = U{t2,tM)U{tM ,ti) where the first 
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factor is anti-time ordered, and the second factor is time 
ordered. Together, symbolicaUy, we can write 



U{t2,Ti) = TcCXp 



Hrdr], {t2^ti), (35) 



where we assume T2 succeeds ti on the contour. Tc is an 
order super-operator that orders the operators according 
to a Hnear order on the contour from earher to later when 
read from right to left. Of course, this makes sense only 
when the exponential function is expanded as a sum of 
polynomials of Hr- The integral is defined as contour 
integral, which we'll describe further in the next section. 
U has a group property on the contour, i.e., if T2 >- 
Ti, then U{t3,T2)U{t2,ti) — U{t3,ti). In addition, if 
Ti ^ T2, we define U{ti,T2) — U{t2,ti)^^ . With the 
evolution operator defined on the contour, we can define 
Heisenberg operator on the contour as 



0(T) = C/(r,4)-iOC/(r,i+). 



(36) 



This definition agrees with the usual Heisenberg operator 
and is independent of the branch index a as an operator 
acting on a vector of Hilbert space, if the Hamiltonian is 
independent of branches, which normally is. However, if 
0{t) is under the contour order sign its position is 
dictated by the contour variable r. 

We are now in a position to define the contour ordered 
Green's function, as 



G(t,t') 



h 



Tr 



p{to)TcU{T)u{T ) 



(37) 



where u(r) now stands for a column vector with Uj{T) as 
elements, and the superscript T stands for matrix trans- 
pose, so that G is a square matrix. 

Working in the branch component form, G(r, r') — 
Gir,t"'') = G""' {t,t'), we obtain four different Green's 
hmctions. We can identify these Green's functions with 
the ones defined earlier by comparing the meaning of 
contour order operator and time-order, anti-time order 
operator. If both rs are on the upper branch, contour 
order is the same as time order, so we have G'^^ = G*. 
Similarly, if both are on the lower branch, contour order is 
equivalent to anti-time order, G = G'. However, if r is 
on the lower branch, and r' is on the upper branch, then 
we don't need to swap positions for all values of t and t' 
for the operator u, so the definition of contour order is 
equivalent to the greater Green's function, G ^ = G^. 
Similarly, G"' = G^ . We can write this as a 2 x 2 matrix 
in the space of branches as 



G 



G+^ 


- G+- ' 




G< \ 


G-^ 


- G— 


Ht 


G* ) 



(38) 



Consider coupled harmonic oscillators (defined by the 
Hamiltonian H , Eq. (|61l) below), initially in some mixed 
state p{to) for which we assume Wick's theorem is valid. 
If the system is then driven by an external force, F{t), 



which is branch dependent, i.e., with an additional r- 
dependent external potential V{t) = —F{t)'^u{t), how 
the density matrix will change? Schwinger gave a result 
[l|, (in our notation) 

Tr{U{tM,to)p{to)U{to,tM)) = (39) 

exp (^-^ J^J^F{TfG{T,T')F{T')dTdT'^ , 

which motivated him to introduce the G°'°' . 



IV. CALCULUS ON CONTOURS; 
CONVOLUTION, TRACE, AND DETERMINANT 

This section is for the mathematically inclined read- 
ers. Those interested in applying the Green's functions 
to physics problems can skip this part in the first read- 
ing. A hallmark of NEGF is the contour valued function. 
To be able to work with the contour functions, we like 
to make a few remarks as how differentiation and inte- 
gration are done on the contour. Analogous to calculus 
on the complex plane the derivative on the contour is 
defined in the usual way. 



hm 

dr At^o 



At) - /(t) 
At 



(40) 



where the function /(r) is equivalent to two functions, 
f^{t) and f~{t), for the upper and lower branch, re- 
spectively. On the upper branch, the definition coincides 
with the usual meaning of derivative with respect to t 
with At = At. On the lower branch, the situation is the 
same, as whether At is positive or negative, it always re- 
sults as the derivative with respect to t. So symbolically, 
we say 



A. 

d7 



d_ 

It' 



df{r) ^ dr{t) 



dr 



dt 



(41) 



The integration is defined very much like contour inte- 
gral. 



dr = 



dt^ 



dt' 



E 



adt 



(42) 



cr = ±l •'lO 



The plus and minus signs on t are to make sure the in- 
tegrand function takes the proper branch indices. If the 
integrand is independent of the branches, then the value 
is zero. 

We define the 9 function on contour as 6{t,t') = 1 if 
T >~ t' where t succeeds t' on the contour, and other- 
wise. We define the 6 function by J(t, t') — d9{T, t') / dr. 
One can convince oneself that 



5{t, t')^ a6a^„,d{t-t'), 



a, <T 



±1, 



(43) 
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where 5 a, a' is the Kronecker delta and 5{t — t') is the 
Dirac delta. The S function on contour has the expected 
property that 



<5(r,r')/(r)dr = /(r'), 



(44) 



if a contour contains the point r', and otherwise. 

In an NEGF calculation, e.g., in collecting terms to 
form a Dyson equation [e.g., Eq. (|8T|) below], one often 
encounters convolution of a certain type on the contour. 
In the theory of full counting statistics, one also needs to 
evaluate trace or determinant defined on contour. Part 
of this is presented in Ref. [2^ in appendix C. In the rest 
of this section, we address these issues, but first some 
notations: 

'A' will mean matrix function with contour times, i.e., 
A — >■ Ajji(T,T'). A{t,t') denotes a matrix with ele- 
ments Ajji with contour time variables explicitly speci- 
fied. A'^'^'{t,t') are the components of A. A''"' = a A""', 
which has the effect of flipping signs for the bottom two 
entries of matrix A, i.e.. 



A 



A* A< 

-A> -A* 



(45) 



A is defined as a 45° rotation in the space of branches 
from A. With the help of Pauli z matrix. 



1 

-1 



and the rotation matrix 



we define 



A^R^ a,AR = R^ AR. 



(46) 



(47) 



(48) 



This is known as a Keldysh rotation (other conventions 
are also used, e.g.. Rammer [35|, Chap. 5.3). For any 
A"'^ , the effect of the Keldysh rotation is to change to 



A = 



A"- A^ 
A^ A'' 

l(A'-A<+A>- A\ A* + A* H 

2[a' + A*-A<-A>, A<- a* 



(49) 



A<+A> 



We should view the above as defining the quantities A^, 
A"-, A^, and A'^ . In particular, A^ ^ A< + A> , as one 
usually might expect, but is equal to {A* + A* + A^ + 
A>)/2. We call A^ the Keldysh component. Although 
A^ need not be 0, it is still true that A'' - A" ^ A> - A< . 
For the Green's functions satisfying the relation ([3T|) we 
get 



G" 



(50) 



The G^ component is due to the relation among 
Green's functions. 

Convolution of n contour matrix objects is defined as 

AB--D = 

J dT2dT3 ■ ■ ■dTnA{Tl,T2)B{T2,T3) ■ ■ • i:)(r„ , T„+l ), (51) 

where the usual matrix multiplication in the indices j is 
implied, and the first and last variables are left free. So 
the result is a matrix function of ti and Tn+i- 

We note that matrix equations are invariant under the 
Keldysh rotation defined by Eq. (|48p . In the normal situ- 
ation when the Green's functions in the Keldysh-rotated 
space are block-upper triangular, the convolution in real 
time or product in frequency domain is still upper trian- 
gular. 



C"" 



A"- A^ 
A" 







(52) 



Multiplying through the matrices, we find C — A^B^ 
and similarly for the advanced component, as well as 
= A^'B^ + A^B". One can also show that G<'> = 
yj^r^<,> _^ j^<,>j^a ^gjjjg ^\yQ general relations among 

the Green's functions. These results are known as Lan- 
greth theorem 38]. Using this technique, it is also fairly 
easy to find the component form of the Dyson equation, 
G = g + gT.G, as 

G"- = / + .g'■S'■G^ (53) 
G'"^ = + /S''G^ + .g'^E^G'^ + g'^T.^G". (54) 

Explicit solutions can be written down as (Ref. Issl . 
Chaps. 4 and 5) 



G'^ = {{g'^y^ - S*-)-!, 



(55) 

G< = (1 + G'^E'').g<(l + S'^G") + G''S]<G''. (56) 

The last equation above is known as the Keldysh equa- 
tion. The first term in Eq. (fSBj) is if {g^)~^g^ = 0. 
This is the case for ballistic systems in steady states. 

Back to the contour ordered version of functions. The 
identity 1 (in the combined contour time and matrix in- 
dex space) is defined by the requirement lA ~ A, or 
AlB = AB. As we can see we can define the iden- 
tity as I5{t^t') where I is the identity matrix while 5 
takes care of the contour space. The inverse of A is de- 
fined by AB = BA = 1, B = A~^ . The inverse has 
the usual meaning if we represent it in the discretized 
version of contour time, then AB = I, B = A~^, 
where the tilded versions are the usual matrices, de- 
fined by discretizing the time with a uniform spacing At 
and indexed by a triplet of cr, j, and ti{— iAt), and 
A-aju-.a-'j't^, = aA'^p {ti,U')At. The tilded version A is 
useful for numerical computation. 

A trace on the contour is defined by integrating over 
all the contour time Ti and the usual matrix trace, which 
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can be represented in a number of equivalent ways: 
= E/ 



= Tr(i) = Tr(^) = Tr(A). (57) 
Finally, the determinant of A is defined through trace, 

„Tr In A 



det(A) = e^''"^ = ^lini^det(A). 



(58) 



With this definition for the determinant, we see that 
det(/(5(r, t')) = 1. If A has the form 1 + M, we can com- 
pute the logarithm of the determinant through a sum of 
traces. 



lndet(l + M) = Trln(l + M) 



Tr ( A/ - -Af2 + 



(59) 



This formula has been used to obtain practical numerical 
method for the computation of cumulants in full counting 
statistics [23. 



V. EQUATION OF MOTION ON CONTOUR 

The equation of motion method is a simple and con- 
venient way to get started in an NEGF calculation. In 
order to deal with the set of Green's functions defined in 
real time, which completely characterizes the system, it is 
useful to consider the equation of motion on the contour 
C [l^. For a ballistic system, such equation closes, so a 
complete, exact solution is possible. However, there are 
some subtleties (on the boundary conditions/initial con- 
ditions) as to how these equations can be solved |;39|] . The 
approach taken here is to express the unknown Green's 
function with what we know, e.g., decoupled equilibrium 
Green's functions discussed in earlier sections. For in- 
teracting systems, the equations become an infinite hier- 
archy, similar to the Bogoliubov-Born-Green-Kirkwood- 
Yvon (BBGKY) type of equations. These hierarchical 
equations can be put into an integral form which is then 
equivalent to the Feynman-diagrammatic expansion of 
the problem. 

The starting point to obtain the equation is the calcu- 
lus rules outlined in Sec. IV, and a generalization of the 
Heisenberg equation of motion on the contour. 



(60) 



where an arbitrary operator is defined on contour ac- 
cording to Eq. (pG]) . Both the derivative and the operator 
with the contour variable r are equivalent to the ordinary 
derivative with respect to time and Heisenberg operator 
at t as far as the effect of operator acting on Hilbert space 



is concerned. So the above equation is totally equivalent 
to the ordinary Heisenberg equation of motion. The only 
difference is that when under the contour order sign, T^, 
the position of the operator needs to be at a proper place 
ordered according to the contour time t. 

We illustrate the idea of the equation of motion method 
with a simple example. Consider a system of coupled 
harmonic oscillators, not necessarily in equilibrium, with 
the usual Hamiltonian 



(61) 



where K is a, symmetric, positive definite spring constant 
matrix, u is a column vector with component Uj, and p 
is the conjugate momentum vector. Since the transfor- 
mation from Schrddinger to Heisenberg operator defined 
on the contour is a unitary transform, the commutation 
relation holds for the contour variables at equal time. 



(62) 



where / is an identity matrix having the size equal to the 
number of degrees of freedom of the problem. 

We write the contour ordered Green's function of the 
full system in terms of the 6 function to facilitate easy 
differentiation, 

G(r, r') - - jTr [p(to)T,«(r)«(T')^] (63) 

= ~^0{t,t') {u{T)u{T'f) - ^9{T\r) {u{T')u{Tff . 

For notational simplicity, we use angular brackets to 
denote average over the density matrix, i.e., (• • •) = 
Tr[p(io) • • •]. We now differentiate with respect to t. 
There are two places that depend on r, one in the 9 
function and another inside the average on u. Using 
de{T,T')/dT = S{t,t'), de{T',T)/dT = -(5(t,t'), the 
Heisenberg equation du{T)/dT = u{t) ~ p{t), and the 
Leibniz rule, we find 

dG{T,T') ^ irr ^ I ,vr\ 



dr 



--<5(r,r')(Kr),u(rr]) 



(64) 



We have combined the two terms proportional to the 9 
functions as a contour ordered one, and combined the two 
terms of derivatives of the 9 functions as a commutator. 
Since (5 is unless t — t' , we can set the second argument 
to T. But equal time coordinates commute, so the second 
term is [Actually, it is oo • 0, but it is safe to set it to 0]. 
For phonons, first order equation does not close, hence 
we take one more derivative to obtain. 



d^G{T,T') 



^-S{T,r'){[u{T)Mr'f])- (65) 

Commuting the Hamiltonian H with u twice, we obtain 
the Heisenberg equation il = —Ku, which has the same 
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form as the classical equation of motion. Further, using 
the canonical commutation relation, Eq. (j62p . the second 
term can be simplified; we obtain the equation of motion 
of a coupled harmonic oscillator system as 



9r2 



+ KG{t,t') = -5{t,t')I. 



(66) 



We consider the application of this equation to the 
problem of thermal transport in a ballistic system. We 
create a nonequilibrium but well-controlled situation by 
partitioning the whole system into three regions, called 
left lead, center region, and right lead. Each one of the 
regions will have a well-defined initial density matrix. 
Thus the matrix K takes the form 



K = 



yCL j^C yCR 

V"^ 



(67) 



where the submatrices if", a = L,C,R are symmetric, 
and V^^ = {V^^ f, V^^ = (t/^^)"^. The sizes of the 
matrices are considered finite at the moment. If we like to 
obtain the steady-state result, we'll send the sizes of the 
leads to infinite at the end of the calculation. In terms 
of the Hamiltonians of the subsystems, we may write. 



H = Hl+Hc + Hr + uiV'-'-uc + ujjV'-uc, (68) 

where the last two terms correspond to the interaction 
of the leads with the center, and the decoupled systems 
have Hamiltonians, Ha = {l/2)p'^pa + {l/2)u^K°'Ua, 
a = L,C, R. 

We split the K matrix into diagonal and off-diagonal 
terms as K = D + V, where 



D = 





K'^ 



lv= V 






yLC 





CL 





yCR 





yRC 






(69) 

With this split of the decoupled ones and interaction, it is 
easy to verify that the following Dyson equation is valid. 



G{T,T')^g{T,T')+ f dT"g{T,T")VG{T",T'), 

Jc 



(70) 



where the contour C is the standard contour from 
to tM and back to . We should view the functions G 
and g strictly defined only in the time interval [to,tM]- 
The advantage of working on this interval instead of the 
Keldysh open domain (— c», +oo) is that we can treat the 
transient as well as steady state on an equal footing. The 
small g given above is defined by 



d'gir, 



9r2 



DgiT,T') = -S{T,T')I. 



(71) 



Symbolically, we write G — g + gVG = g + GVg, where 
the multiplication should be understood as a convolu- 
tion on the contour. Equation (j66|) is obtained if we act 
the differential operator, Id^/dr^ + D, on both sides of 



Eq. (|70|) . using Eq. (|7T|) and the property of the S func- 
tion. The Dyson equation fulfills our goal of expressing 
the unknown, possibly nonequilibrium Green's function 
G in terms of simpler Green's function g. 

However, it is not a good idea to focus on solving the 
differential equation (|7ip as the solution is not unique. 
If g' satisfies g' + Dg' = 0, then g + g' also satisfies 
Eq. (|7ip . Thus, we'll have to fix the small g accord- 
ing to their original definition, Eq. (j37p . using the ini- 
tial density matrices, with the decoupled Hamiltonian 
h = {l/2)p^p+{l/2)u^Du = Hl+Hc + Hr. Since the 
Hamiltonian is quadratic, with a product of equilibrium 
initial states, the three systems are completely decoupled, 
we have for g at aa' subblock. 



5aQ'(T, r') = --TT[pL{ta)pc{tQ)pR{to)TcUa{T)Ua'iT') 

n 



6a,a'ga{T,T'), a,a'^L,C,R, 



(72) 



where we have {ui^,uc,U]i)'^ — u, and the time evo- 
lution is according to h. Obviously, Eq. ([7T|) is sat- 
isfied with this definition. If Pa{to)s are the equilib- 
rium distributions, then we have established our goal 
of 'building' the nonequilibrium Green's functions from 
initially known equilibrium ones. This will be the 
usual case, but if the initial state p(io) is arbitrary, 

t" —t^ 

an extra surface term, g{T,T")dG{T" ,t') / dr"] ° — 

T — 

%(T,T")/ar"G(r",r')rIl*°+, needs to be added to the 
right-hand side of Eq. (170]) . 



VI. FEYNMAN DIAGRAMMATICS 

The Feynman-diagrammatic perturbation theories are 
the standard techniques to treat interactions in a sys- 
tematic way. We'll refer to the literature for details 
[H m m, ii, IMIil • Most of the earlier literature treat 
systems at absolute zero temperature (e.g., quantum field 
theories). For finite temperature in thermal equilibrium, 
it is the Matsubara formalism that is employed. For- 
tunately, the diagrammatic structures are all the same, 
whether it is nonequilibrium contour order, or T = 
time order, or Matsubara order. 

As an illustration, we study the problem of a ballistic 
system divided into three regions with the Hamiltonian 
given by Eq. (j68p with one additional term for a quartic 
nonlinear interaction which appears only in the center. 



4 ^ 

ijkl 



(73) 



An important step for a perturbative expansion is to sep- 
arate the system into a solvable one and a perturbation. 
A two-step adiabatic switch-on may be used as illustrated 
in Ref. 15. Here, we'll consider a sudden switch-on at 
time to with the decoupled system h = Hj^ + He + Hj^ as 
the unperturbed one and the lead-center couplings and 
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the nonlinear interactions as perturbation. In the interac- 
tion picture with respect to h, the operators and density 
matrix are transformed unitarily from the Schrodinger 
picture by 



(a) 



Oi(t) = e^('-*°)''Oe-^(*-*«)^, 

//,(^^i(«")+^^I"(*"))dt" 



(74) 
(75) 

t > t'. (76) 



So the operators foUow a 'free' evolution of the noninter- 
acting system, and the density matrix evolves through 
the evolution operator 5* according to only the interac- 
tion part of the Hamiltonian. Both O and S can be 
generalized to be defined on the contour (by saying that 
the time t has an additional contour branch index a). 

We would like to compute two quantities. The first is 
the "partition" function or generating function 



= Tr 



(77) 



where the second line is in the interaction picture with 
Viir) = u£(T)^y^^u^(r) + u^^iT)^V^^u^ci^) due to 
the interaction picture transformation. This quantity Z 
is clearly 1, by definition. But we'll take the point of 
view that the full Hamiltonian may be contour depen- 
dent, then Z = Tr[p(to)[/(<Q , tj)] may not be 1, a useful 
point of view when we discuss full counting statistics. In 
addition, we consider the functional form of Z in terms of 
the Green's functions g and look for the relationship be- 
tween the full Green's functions G and Z. The diagrams 
generated in Z are known as vacuum diagrams as there 
are no external lines. Of course, our main focus is the 
second quantity, the contour ordered Green's function, 
Eq. ((37|). When transformed into the interaction picture, 
we have 



Gaa' (t, t') 



-Tr 



pito)Mui{T)ui,ir'f 



(78) 



Jciyiir")+H'^{r"))dr'^ 



— , a^L,C,R. 



There are a number of well-known facts or theorems 
which will be helpful in the development of the Feynman- 
diagrammatic expansion. We'll list them here without 
proofs. 

1) The Wick theorem. This theorem enables one to 
express the product of terms when the exponential is ex- 
panded, in terms of simpler known Green's functions, 
e-g-, 

(reu(l)u(2)u(3)u(4)) = (r,w(l)u(2)) (Tcu{3)u{4)) 
+ {T,u{l)u{'S)) {T,u{2)u{4)) 
+ (r,u(l)u(4))(T,?/(2)M(3)), (79) 

where for notational simplicity, we have lumped a set of 
indices and contour time argument as one single number. 



L 



O -O 

R 



In Z= J + - <0, + 



(c) 



(d) 



(e) 



Gcc — 



lO -f} - 



+ 
+ 



Gcc = = +3 



lnZ = 



+ 6 




FIG. 2. Feynman diagrams for the nonequilibrium transport 
problem with quartic nonlinearity. (a) Building blocks of the 
diagrams. The solid line is for gc, wavy line for g^, and dash 
line for gn; (b) first few diagrams for InZ; (c) Dyson series for 
the ballistic system Green's function Gq^; (d) Full Green's 
function Gcc , and (e) resumed In Z where the ballistic result 
is InZo — — iTrln(l — gcS). The number in front of the 
diagrams represents extra combinatorial factor. 



e.g., (1) = {ai,ji,Ti). The validity of the Wick theorem 
relies on the fact that Inp(io) is a quadratic form in the 
dynamic variables. Each graph comes with a numerical 
prefactor, which can be found by working out a combina- 
torial problem of how many ways one can get a topolog- 
ically equivalent graph due to the Wick decomposition. 

To work out the diagrams, for our problem of the 
center-lead couplings and the quartic nonlinear interac- 
tion, we have several building blocks. First, the pairs of 
u give the decoupled Green's functions ga, a = L,G,R. 
This will be drawn as wiggle, straight, and dotted lines, 
(see Fig. HJa)). These lines are connected possibly in 
two ways, by the V^'~' or V^'~' vertices, which connect 
the center line with the lead lines, and by Tijki which 
connects four center lines. For G we have two external 
terminals labeled 1 and 2. For Z all variables are dummy, 
and need to be summed. 

2) Cluster decomposition theorem, factor theorem. The 
graphs of Z contain connected and disconnected pieces. 
The cluster decomposition theorem is a very general the- 
orem which says that if we take the logarithm, then In Z 
is given only by the connected graphs (see Fig. [^Jb)). A 
similar statement holds in the Mayer's cluster expansion 
in equilibrium statistical mechanics for interacting gases 



10 



(Friedman [4J|, Chap. 6). In addition, the disconnected 
pieces do not enter into the diagrams for the Green's 
functions G. This can be understood in two ways, first is 
that due to the denominator Z in Eq. (|78l) . the discon- 
nected pieces (which contains only vacuum diagrams) get 
exactly cancelled. Alternatively, if Eq. ((3T|) holds, then 
all the vacuum diagrams are numerically 0, so such dia- 
grams do not appear and Z = \ [45]. 

3) The Dyson equation. Certain diagrams can be re- 
grouped and easily summed. Let us first consider the case 
where the nonlinear interaction vanishes. Then the dia- 
grams for the corresponding Green's function can 
only be a linear chain, with binomial combinatorial ways 
of putting the left or right lead lines, see Fig. Efc). If we 
define the self-energies as 

I]„(t,t') = l^^"<?a(T,T')F"'^, a = L,R, (80) 

the terms can be expressed in a recursive way. We thus 
have (for the center part of the full ballistic Green's func- 
tion) 

Gcc = 5c + gc^gc + gc^gc^gc h — 

^gc + gc^G^c. (81) 

where E = + Sj^ and convolution on the contour 
is implied. This is the Dyson equation for the central 
region, which has a similar mathematical structure as 
Eq. (EOD. 

The nonlinear part is more complex, but in terms of 
Gqq (double line in the graphs), the diagrams for InZ 
and full nonlinear Green's function Gcc can be simpli- 
fied, as shown in Fig. [2Ie,d). Finally, we can define the 
nonlinear self-energy as part of the diagrams where it is 
not singly connected (double or more connectivity) and 
with the two external legs chopped, thus giving 



cc 



GO 



cc 



Gcc^nG, 



cc- 



(82) 



4) One can introduce a vertex function (and Hedin-like 
equation) and encapsulate the diagrams more compactly 
using functional derivatives [s^ Hy, \^ . How useful they 
are for practical calculation remains to be seen. 

5) Connection between vacuum diagrams and Green's 
function. This fact seems less well-known. We notice that 
the vacuum diagrams given in In Z and the graphs formed 
by Tr(GccS) are the same, where the trace means both 
for the space index j and contour time r, and the expres- 
sion means to close the two external lines with one more 
self-energy line. However, the combinatorial prefactors 
differ. This difference can be removed if we differenti- 
ate with respect to the self-energy E. Thus we have the 
following identity. 



SlnZ = -Tr (GccS^) : 



(83) 



where the variation S means the functional form of the 
self-energy is varied while gc holds constant. This re- 
lation can be derived in a more rigorous, algebraic way 



VII. LANDAUER FORMULA 

So far we have studied the properties of the Green's 
functions and how such functions can be calculated for 
general linear or nonlinear systems. In this section, we 
look at one of the most important physical observables 
in transport, i.e., the thermal or energy current. The 
energy current transported out of the left lead is defined 
as 

'dHLit)' 



hit) 



dt 



(tti(t)^F^^uc(t)), (84) 



where the angular brackets denote trace over the initial 
density matrix p{to) and the operators are in Heisenberg 
picture at time t. This energy (per second) is presumably 
transferred to the center or the coupling between the left 
lead and the center, since the energy of the whole system 
is conserved and the left lead is connected directly to the 
center but not to the right lead. 

We need to connect the definition of the current with 
the Green's functions. Using the definition of G^^ or 
Gqj^ (or G^^/2) we can write 

d 



IUt)=^h—^r[G<fit,t')V^^] 



(85) 



The trace above is in the sense of an ordinary ma- 
trix trace by summing over the diagonal elements. The 
Green's functions above use the mixed lead and center 
degrees of freedom. Observing the fact that centers are 
usually more complex but finite, and leads are simple 
(free phonons) but may be infinite, we can try to relate 
Gcc to GcL using the Dyson equation, ([70)1 . in the form 
G = g + GVg. Working out the GL component of this 
block matrix equation, we obtain Gcl — GccV^^ gL, or 
in full detail 



Gcl{t,t')^ / Gccir,T")V^'^gLiT",T')dT". (86) 
Jc 

Although we obtained this equation from the ballistic 
Dyson equation, the specific properties of the center has 
not been used. It turns out that this equation is also 
valid if the center is nonlinear. This can be shown by 
looking at the equation of the motion of Gcl directly. 
Using the Langreth theorem for the Gqj^ component and 



substituting the result into Eq. 



CL 

we obtain 



lL{t) = ihTr 



G^cc{tX)g^^tit",t')\ 



G<^{tX)g^mt",t')\ 



dt", (87) 



where we have used the definition of lead self-energy, 
Eq. (HOI). 

The above formula is valid for any time t. If steady 
state result is required, we can send to — — oo, and per- 
form a Fourier transform and after applying the convo- 
lution theorem for Fourier transform, we obtain 



— — TiwTr 



G^[lo]^<[c.] + G<[lo]^1[co] 
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n+oo_ 



An' 



hoj Tr 



G<S> 



G>I]< 



For simplicity, we have dropped the CC subscript on G; 
we have also omitted the w arguments in the second line. 
The last equation can be obtained by taking (II + /£)/2 
since II must be real, and using the relations among 
Green's functions discussed in Sec.lIIII The above formula 
is known as the Meir-Wingreen formula first derived for 
electronic transport (lo| . 

The Meir-Wingreen formula is valid for ballistic as well 
as interacting centers. If the center is ballistic, i.e., I]„ = 
in Eq. the formula can be further simplified. The 
result for the ballistic system is called Landauer formula 
with a transmission function T{uj) = Tr(G''riG"rfl) 
known as the Caroli formula 0. Thus, using the defi- 
nition for the lead spectral function 



(89) 



and the relations (as a consequence of the Dyson equation 
for the center, see, e.g., Datta [4§|, Chap. 8) 



G< = G''(E< E<)G^ 

we finally obtain 



It. 



2^ 



(90) 
(91) 



(92) 



where fa = l/(e'^°''" — 1)^ a = L, R, is the Bose-Einstein 
distribution for the leads. A very detailed derivation of 
the above is given in Leek [47'|, Chap. 3. A similar for- 
mula for electrons was first given by Landauer [s^, [Hi 
from a wave scattering point of view. It is appropriate to 
call Eq. (j92|) Landauer-like formula and it has been de- 
rived in a number of different ways for thermal transport 

As it is seen that the surface Green's functions and 
gu are the important inputs for the noncquilibrium trans- 
port problems, it is required to develop algorithms to cal- 
culate them efficiently for realistic systems. Fortunately, 
this has already been done quite early [6l|, HI] for elec- 
tron transport, but it applies equally well for phononic 
systems. Algorithmic procedure for the calculation of 
the surface Green's function and thus the self-energies 
are reviewed in Ref. [isl 

In the rest of this section, we give a simple example of 
calculating the Green's functions and transmission coef- 
ficient for a uniform one-dimensional chain with the force 
constant matrix K which is 2k + ko along the diagonal 
and —k along the two first off-diagonals. We split the 
system into three regions with Nl, Nq, and Nn number 
of particles for each. The first step in such a calculation 
is to determine the surface Green's functions. The eigen- 
values and eigenvectors of the uniform tridiagonal matrix 
can be obtained analytically [H, , given 



2k 1 



cos 



(qn)) + ko, 



N+1 



(93) 



sin(g„j), n = l,2,---,iV, (94) 



where N can be one of the Na, a = L,C,R. 
We construct an orthogonal matrix S, S'^S — 
I, hy S ^ {u\u^,---,u^) such that S'^KS = 
diag(r2^, r22, • • • , ^%)- Each mode follows the results of 
the single degree harmonic oscillator discussed in Sec. II. 
Then the retarded Green's function in the time domain 
for a chain of N sites is 



g-{t)=Sdiag\-e{t) 



sin(57ji) 



(95) 



For steady state, we need to have an infinite lead, N — 
oo. In order to take this limit, it is more convenient 
to solve the retarded Green's function in the frequency 
domain, 

((w + zry)2-if"),g^[w] =/, a^L,C,R. (96) 
The (1,1) element can be found analytically, given 

A 1-A2^ 
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(97) 



k 1 - \™+'^ ' 
where A is the solution of the quadratic equation, 

k\-^ + {u + ir]f -2k-kQ + k\^{), (98) 

with |A| < 1 (the small imaginary number irj makes this 
choice unambiguous). Since |A| < 1, in the limit N ^ oo, 
we obtain a simple result for the surface Green's function 
of the semi-infinite lead as g^[u)]ii = —X/k. Using the 
definition of self-energy, Eq. (1501) . we obtain for the ma- 
trix elements = Y,'^f;[uj]N^N^ = —kX, and for 
all other elements. 

The retarded Green's function of the coupled system 
in the center can be obtained by solving the Dyson equa- 
tion, Eq. (|8T|) . of the retarded component, using Langreth 
theorem and in the frequency domain. 



(99) 



However, we can also obtain G'' by considering a similar 
equation as for g for the whole space domain, with integer 
index j, I vary from — oo to -|-oo, i.e., 

fcG^^_i,, + ((c^+ir;)2-2fc-fco)G,-,-f fcG^Vi,/ = (100) 

Since G'' must be translationally invariant in space in- 
dices, and must decay to when j or Z — ±oo, we have 
[H Gji[uj] = cAl^-'l, where c = 1/[(A - 1/A)fc] is fixed 
by the j = I diagonal equation. Finally, the transmission 
coefficient is found by the Caroli formula. After some 
algebra, one finds T{uj) 



1 if vAJq < uj < ViF+fco and 
otherwise. Of course, this simple result is expected if 
one thinks of it from a wave scattering picture (66l . \qt\ 
without doing any calculation. 
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VIII. MULTIPLE LEADS, LEAD-LEAD 
INTERACTION 

In this section, we give some formulas without much 
derivation for balhstic transport in some more settings. 
The first is that of multiple leads. This is much in parallel 
to Biittiker's theory [i^, l68l - [7l| for electron transport 
with multiple leads. We define 



(101) 



for the transmission coefficient from the a to a' lead. 
Then the current out of the a lead is given by summing 
over contributions from all the other leads, 



Ia= f '^Tt^ y2 T'a'ai^)ifa- fa')- 

Jo '^^ lU- 



(102) 



Three-terminal problems are studied in Ref. [T^ and [TJ 
where one of the terminals is treated as a Biittiker probe, 
i.e., the third terminal is required to have zero current, 
determined self-consistently by adjusting its bath tem- 
perature. This mimics an inelastic scattering, thus result- 
ing in thermal rectification even for ballistic systems, as 
well as diffusive transport for long chains |74l - l76l | (where 
each atom gets a probe or self-consistent reservoir). How- 
ever, the effective decoherence is a bit artificial, and 
its relevance to truly nonlinear systems is not clear. A 
four-terminal problem is treated in Ref. [t^ for the spin- 
phonon or phonon Hall interaction (still a ballistic prob- 
lem since the "interaction" is bilinear in coordinates and 
momenta). 

Now we come back to the two-lead problem again. In 
the standard modeling of such systems, one always as- 
sumes that there is no interaction between the left lead 
and the right lead. This, of course, can be achieved if 
the center region is large enough and the interactions 
are short-ranged. However, in practical calculation, one 
always finds a small residue of the left-right lead interac- 
tion. Now the question is, can we have a generalization 
of the Caroli formula so that the left-right lead interac- 
tion is allowed? The answer turns out yes, with a new 
formula [78]: 

T,(c.)-Tr(G?,ifiGl«f«), (103) 
where the new tilded lead spectral function is 



Ta = i 



a^L,R. (104) 



The subscript red indicates that the surface Green's func- 
tions are not the full one but reduced subblock large 
enough so that outside that sizes, the left-right couplings 
(as well as center-left, center-right) are 0. Also, the ma- 
trix G£^ is not infinite large but a finite piece consistent 
with the sizes of (/red, a- The retarded Green's function 
in the LR subblock can be obtained by solving a Dyson 
equation 
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where V^^ = (V^^)'^ is the left-right lead coupling ma- 
trix and G^2. — (^lrV ■ is also possible to give an 
expression as a special case of the above for the trans- 
mission where the system has no center and left and right 
leads are directly connected (an interface problem), as 



(107) 



where jl = iV^^{g}^ - g1)V^^. The starting point of 
deriving these results is Eq. ([TO)) , where V is non-zero 
for all the subblocks except on the diagonal. We refer to 
Ref. (78) for more details. 



IX. FULL COUNTING STATISTICS 

In thermal transport, the current is the first and most 
basic quantity to look at. However, other related quan- 
tities are also relevant and important. One of them is 
the current fiuctuation. For example, we may consider 
the time-displaced current-current correlation function 
(much work has been done for electron shot noise [tH). 
In a transient situation, the fluctuations are large. Due to 
the stochastic nature of the baths, each individual exper- 
iment will give a different result Q (— Ql) if we measure 
the amount of energy transferred in a fixed amount of 
time tjv/ out of the left lead. Thus, it is interesting and 
useful to look at the distribution of the energies. This 
distribution satisfies certain 'fiuctuation theorem' under 
certain conditions (e.g^ long time), which is now a very 
hot area of research [79l - [8l| . The complete distribution, 
or equivalently the moment generating function (charac- 
teristic function), 



z{0 



+ 00 



dQe^««P(Q), 



(108) 



reveals more about the system, especially its quantum 
nature. When the function Z{£) is known, the moments 
of Q can be computed by taking derivatives, (Q") = 
(?"Z/9(z^)", and then setting C = 0, and the cumulants 
are defined through ((Q")) = 5" In Z/9(iO"l?=o- In par- 
ticular, the first moment or cumulant is proportional to 
the current in long time, {Q) = {{Q)) ~ tul] the second 
cumulant is the variance {{Q'^)) = (Q^) — (Q)^, and so 
on. This problem of the study of the distribution of the 
transferred quantity is known as full counting statistics in 
the electronic transport literature. There, the number of 
electrons transferred in a given time is a discrete quan- 
tity, thus the word 'counting' is appropriate. Phonons 
cannot be counted (however, see Ref. |82 ). and also we 
are not interested in the number of phonons since it is 
not a conserved quantity. What we do here is to measure 
the amount of energy, a continuous variable, transferred 



13 



from the left lead into the center. In statistical mechanics 
literature, Z is related to so-called large deviation prob- 
lem (when time tu approaches infinity) [ssj . 

In defining the generating function, the most impor- 
tant observation is that the energy transferred, Q, is not 
associated with the eigenvalues of a quantum-mechanical 
operator. Instead, it is computed by the difference of the 
energies of the left lead at two different times, Q = a — h, 
where a and b are the eigenvalues of at time and 
tM, respectively. Using this two-time measurement pro- 
tocol and the standard von Neumann's interpretation of 
quantum measurement, we can derive a very general for- 
mula for Z(^) using product initial state given by [tqHs^ I. 



a,b 

= Tr 



(109) 



[pL{to)Pcito)PR{to) 



where P{b, a) is the joint distribution for the event his- 
tory which at the initial measurement resulted in the left 
lead with energy a and b at the second measurement at 
time tM, assuming discrete energy spectrum. If the initial 
state is not a product of equilibrium states, the situation 
is more complicated. We refer to Ref. [25l for more details. 

We can write Z in terms of the modified evolution 
operator Ux{t,t') governed by a modified Hamiltonian 
Hx = e'^-^^i/e"'^-^^ which is contour branch depen- 
dent with X — on the upper (forward) branch 
and a; = on the lower (return) branch, giving Z = 
{U^/2ito,tM)U^^/2itM,to)). Transforming into the in- 
teraction picture with respect to the decoupled system 
h = Hl + He + Hfi, we obtain 



^(0 



(110) 



where = u£ V'^^uc+u^V^^uji in the ballistic case. 
The effect of measurement is to replace the variable ul 
by a transformed one. In the interaction picture, this is 
equivalent to shifting time argument, i.e., u^{t) = u(t + 
7ix(r)) where the amount of shift depends on the branch. 
It is on the upper branch and Ti£^/2 on the lower 

branch. We then use the Feynman-diagrammatic method 
to expand the exponential and group the various terms. 
After some simplification, we obtain [2^, [sH] 



InZ(e) = --Tr,-,ln(l - G^c^t), (111) 

where the trace is over the ordinary space index j as well 
as over the contour time r, G^C standard ballistic 

contour ordered Green's function of the center, while the 
important new self-energy 

S^^(t, t') = Si {T + nx{T), t' + nx{T')) - El(t, t') (112) 

is the difference between the left-lead argument-shifted 
self-energy and the original standard lead self-energy. An 
alternative expression valid also for interacting systems 



[86l . ISTj l for the derivative of In Z with respect to is 
[using Eq. 
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dr / dr'Tr 
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GcciT,T') 



a(*0 



(113) 

We can think of this equation as a generalization of 
the Meir-Wingreen current formula to the full counting 
statistics. The meaning of the tildes there is that we ex- 
press all quantities in terms of small ga, a — L,C, R, and 
then replace all occurrence of by 



~9Lir,r') = 



(114) 



which is simply an argument-shifted version of the origi- 
nal left lead Green's function. 

We note that Eq. (|llip or (|113p is defined on the seg- 
ment of the contour C, so the result is valid both for 
transient and steady state (if we take to — > — oo and 
tM ^ +oo). In the long-time limit, after transforming 
the Green's functions into frequency domain using the 
property of time-translational invariance, and employing 
the standard relations of Green's function and some al- 
gebra, the cumulant generating function for a ballistic, 
left-center-right junction system is then given as 



InZ(e) 



-tM 



+ 00 



duj 
An 



/l(1 + /fl)(e'«"" - 1) + Ml + /L)(e-'«''" - 1)1 .(115) 



This is the phonon analog of the famous Levitov-Lesovik 
formula for electrons [H, H^, first given by Saito and 
Dhar The long-time generating function satisfies 

the relation Z(^) ~ Z(^—£_ + i{f}R — /?_l)), which is a 
form of the Gallavotti-Cohen symmetry [91] . Classical 
versions are given in Refs. 92 and 93. Similar generating 
functions have also been obtained for systems with driven 
forces [15, with the left-right interaction term, uj^V^^UR 
[S^ , as well as an extension to nonlinear systems starting 
from Eq. ((83)) through a counting field ^-dependent non- 
linear self-energy S„ [s^l- Further details can be found 
in Ref. 45. 

We present a numerical example applying the theory 
outlined above for a uniform one-dimensional chain. In 
Fig. [3] we plot the first four cumulants of energy trans- 
ferred for ballistic system with one atom at the center 
which is connected with two finite-size leads. We use 
the formula given in Eq. (Illip and numerically evalu- 
ate the derivatives of InZ(^) with respect to the count- 
ing parameter The self-energy for the finite lead 
is calculated using the lesser version of Eq. ([^ . Gf^f^ 
is obtained by numerically solving the Dyson equation 
given in Eq. (|8ip . The plot shows that for finite leads all 
the cumulants reaches a quasi-steady state with a finite 
recurrence time tr which depends on the length of the 
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FIG. 3. (Color online) Plot of the cumulants of heat {(Q">> 
for n = 1, 2, 3, and 4 for one-atom center connected with two 
finite leads (one-dimensional chain) as a function of measure- 
ment time tu- The black (solid) and red (dashed) line cor- 
respond to Nl = = 20 and Nl = = 30, respectively. 
The initial temperatures of the left, center, and the right 
parts are 310 K, 360 K, and 290 K, respectively. "We choose 
k = leV/(uA^) and fco = 0.1 eV/(uA^) for all particles. 



full system and the velocity of the phonon waves. Af- 
ter tr phonon waves which are scattered back from the 
boundaries interfere and this results in the cumulants to 
oscillate rapidly. Similar results are obtained for a left- 
right lead problem without center [6^ . For infinite size 
leads (94i] complete irreversible behavior emerges and the 
system achieves a unique steady state with infinite recur- 
rence time. The slopes in the quasi-steady state regime 
match with the predicted values obtained from Eq. ()115|) . 



one-point Green's function to the usual two-point one. To 
lowest order of the T^fe cubic nonlinearity, there is only 
one diagram (Feynman diagrams up to second order in 
h are given in Ref . US) • We call it the lollipop diagram 
and is given, algebraically, 

G'.W = Y.^imnJdr'GUr\r')G'njir',r), (117) 



where the superscript refers to the unperturbed, bal- 
listic system Green's functions. In equilibrium or steady 
state, the contour is from —oo to -l-oo and back to — oo. 
Thus, the result is independent of time and branch index. 
The Green's function can be further expressed in real 
time or in frequency domain using the Langreth rules 
(e.g., G^(0) / dtG^{t)). The displacement relative to 
the ballistic equilibrium (which gives (uj) — 0) is then 
computed from iTiGj. Numerical results for nanotubes 
and graphene sheets with Brenner potential are given in 
Ref. 123- It is interesting to note that the thermal ex- 
pansion coefficient in the radial direction of nanotubes 
and graphene sheets are negative at room temperature 
or below. The same method is applied to the study of 
thermal contraction in silicon nanowires f96l], as well as 
multi-layered graphene [97| . 

Another simple application of the lowest order pertur- 
bation expansion of the nonlinear diagrams is the phonon 
life-times. If we work in the normal-mode representa- 
tion so that in the noninteracting system each vibra- 
tional mode is diagonal with the retarded Green's func- 
tion given by l/((w + i-qY ~ ^Iq), then the effect of non- 
linearity can be interpreted as to shift the vibrational 
frequencies of each mode q and to give a damping or 
finite life-time of the mode: 



X. NONLINEAR SYSTEMS, PERTURBATION - 
THERMAL EXPANSION, PHONON LIFE TIME 

NEGF offers a straightforward treatment for a per- 
turbative expansion result when systems are nonlinear. 
We illustrate this with two applications: the problem 
of thermal expansion and phonon life time. The ther- 
mal expansion in bulk systems is usually treated with 
the standard Griineisen theory where the size-dependent 
vibrational frequencies (through the Griineisen parame- 
ter) are the key parameters, see, Ashcroft and Mermin 
[95l |. Ghap. 25. For a finite system, we look directly at 
the equilibrium displacements with a proper boundary 
condition (e.g., one side of a graphene sheet is fixed) in 
comparison with a corresponding ballistic one due to the 
lowest order nonlinear effect [IJj . To this end, it is useful 
to introduce a one-point contour ordered Green's func- 
tion. 



G,(r) 



-(TcUj(t)). 



(116) 



The contour order does not play any role here, but it is 
convenient and uniform in notation when we relate the 
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We'll assume that the effect of nonlinearity is small 
and the retarded Green's function is essentially peaked 
around Og with a small shift by a complex number 
Ag — i/Tg where Tq is the phonon life-time of mode q. 
Under such approximation, we can identify the real and 
imaginary part of the retarded nonlinear self-energy as 

[iiii 



ReS];^,[0g]«20gAg, 



ImS;,g[f^g] 



20„ 



(119) 
(120) 



It is interesting to note that at the lowest order of ap- 
proximation from NEGF, the result agrees exactly with 
the Fermi- Golden rule result. 

We can use the kinetic theory formula to compute the 
thermal conductivity, k — ^CqV^Tq, once the phonon 
life-time is known, where Cq is heat capacity per unit 
volume of mode q, and Vg is the group velocity of mode 
q. However, such approaches are not very rigorous, as 
many assumptions have gone into it. 
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XI. NONLINEAR SYSTEMS, MEAN-FIELD 
APPROXIMATION 



Nonlinear problems are the heart of matter and the 
holy grail of thermal transport. In principle, NEGF 
solved the problem formally by giving the nonlinear self- 
energy S]„ together with the Meir-Wingreen formula for 
current. However, a practical calculation with good accu- 
racy is immensely difficult, particularly if the sizes of the 
systems are large. As a start, we can use lowest order 
perturbative expressions for the nonlinear self-energies. 
It works to some extent for weak nonlinearity and small 
sizes. Our experience seems to indicate that, with just 
the perturbative terms, it is not possible to produce 
correct diffusive transport for large sizes [13, HI]. The 
nex t step is to use self-consistent Born approximation 
[gsI llOCt llOlj l , keeping only the lowest order self-energy 
diagrams (see Fig. 5 in Ref. KB for the self-energy dia- 
grams). Self-consistency means that the ballistic Green's 
function G" is replaced by the full Green's function G. 
Such an approach gives only qualitatively correct results 
(such as diffusive behavior for large sizes, i.e., the current 
decreases with sizes as l/L). One technical difhculty 
in the self-consistency procedure is that the iterations 
may not converge, as the Green's functions are oscillatory 
functions, nonsmooth and ill-behaved. Some approxima- 
tions do not cons erve ener gy exactly, i.e., Il + Ir 0, 
in steady states (l02l Il03j |. Thus, an accurate, quan- 
titatively correct theory for nonlinear quantum thermal 
transport in a large parameter r egio n (system sizes, non- 
linear strength) is still lacking [l04{ . Many calcula tions 
are still based on solving the Boltzmann equations (l05j | 
which is semi-classical (because simultaneous position 
and lattice momentum of phonon distribution are used 
and certain coherent wave nature is neglected [47j). 

Surprisingly, the mean-field theory under certain re- 
stricted condition (quartic nonlinear, small number of 
degrees of freedom) gives very accurate results in com- 
parison with other methods, in partic ular , in comparison 
with the quantum master equation [l06]. The quartic 
nonlinear model, with the potential of Eq. (|73p. is a bet- 
ter model to study as it is stable with proper choice of the 
coefficient Tijki, while the cubic nonlinear term is unsta- 
ble for sufficiently large displacement and always needs a 
quartic term to stabilize the system, although for almost 
all practical systems, the cubic term should be present. 

Our motivation here is to derive reasonably accurate 
equations for the Green's functions without any ad hoc 
approximation for the nonlinear self-energy, but rather 
look at the Green's functions directly. Since the Green's 
functions will form a hierarchy, we need to introduce a 
general n-point Green's function 



G(l,2,...,n) 



-(T,M(l)w(2)---M(n)), 



(121) 



where the number denotes the complete set of space in- 
dex and contour time variables, e.g., 1 means (ii,Ti). 
By definition, G is completely symmetric with respect 



to the permutation of the arguments. For the quartic 
potential. Green's functions with an odd number of dis- 
placement fields is 0, so we only need to consider n even. 
The lowest one is the two-point Green's function. Apply- 
ing the equation of motion method, also taking care that 
we are only interested in the Green's functions involving 
the center degrees of freedom, we obtain 



-(5(ri,T2)(5jij2 



^ X/ '^h33j4j£,^j3jiro32i'''ii'''i^'''ij'''2), (122) 

where / is the identity matrix, is the force constant 
matrix for the center, and E is the self-energy of the leads. 
This equation is exact, and is the first of the BBGKY hi- 
erarchy relating G(l,2) to G(l,2,3,4). In contour time, 
the action of matrix S on G is a convolution, thus 

(SG)(t,t')-/ ^iT,T")G{T",T')dT". (123) 

Jc 

Equation (|122l) can be "solved" or put into an integral 
form, to give 

G(l,2) = G"(1,2)+ (124) 
d3 dA d5 d6 G°(l, 3)T(3, 4, 5, 6)G(4, 5, 6, 2), 

where the contour-time dependent coupling is defined by 
inserting three contour S functions so that all the times 
are synchronized, and G*^ is the Green's function for the 
ballistic system (when the nonlinear term is 0). We can 
carry on to derive equations for G(l,2,3,4), which will 
then involve a 6-point Green's function. At some point, 
we have to close the equations by certain approximation. 
As there is no particular good reason to prefer one ap- 
proximation over the other, such an approach is overly 
complicated and seems at a loss. Perhaps the simplest 
one among all is to stop as early as possible, thus we 
consider 



-)G(1,2,3,4)«G(1,2)G(3,4) 



(125) 

G(1,3)G(2,4) + G(1,4)G(2,3). 



This equation would be exact if Wick's theorem is valid. 
This approximation is amount to the assumption that 
high-order correlations (4-th order) are small, thus tak- 
ing only 2-point correlations should be already a good 
approximation. We like to point out that this is not a 
weak nonlinear approximation, as it is not obtained by 
truncating a perturbation series. The validity of such an 
approximation can only be tested numerically. Putting 
this approximation for the 4-point Green's function back 
into Eq. (jl22l) or (|124p we see that the result is equivalent 
to having a self-consistent nonlinear self-energy, taking 
only the lowest order diagram 

m6{T,T')J2T,fkiGki{T,T). (126) 

kl 
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FIG. 4. (Color online) Comparison between SCMF (solid 
lines) and master equation (crosses) for the one (left) and 
two (right) particle quartic nonlinear model. For the one 
particle case Q.^ = 60.321 meV/{A^u) and Tun = 0.241 
(black curve), 1.2 (red curve), 2.4 (blue curve) [eV/(A*u^)]. 
In case of two particles, Kn — K22 = 60.321, K12 = K21 = 
-30.165 [meV/(A2u)]; Tnn = r2222 = 0.483, T^i, 1,1.2} = 
-^{1, 1,2,2} = -0.241 (black curve); Tun = r2222 = 2.4, 
^{1,1,1,2} = -^{1, 1,2,2} = -1-2 (red curve); Tun = r2222 = 
4.8, r{i,i,i,2} = -^{1,1,2,2} = -2.4 (blue curve) [eV/(A*u'')] 
and the curly-brackets in subscripts indicate all possible per- 
mutations of the indices. The retarded self-energy of the 
leads E;[cj] = ip J„(a;')/(a; - uj')duj' - iJc(uj), with 

Jc{oj) = e^a;/(l -f- uj'^/loI,), L,R,e = 6.0321 meV/(A2u), 
hujD = 10 eV, Tl = 1.25r and Tr = 0.75T, which corre- 
sponds to the Lorentz-Drude model of heat baths. 



dure is probably one of the key ingredients for treating 
strongly anharmonic systems within the NEGF frame- 
work. 



XII. REDUCED DENSITY MATRIX FOR 
BALLISTIC SYSTEMS 

For ballistic systems, we can calculate the n-point 
Green's functions (for the center degrees of freedom) for 
any values of n. Clearly, the set of all Green's functions 
completely characterizes the steady state of a nonequi- 
librium system. In fact, only two-point Green's function 
is needed as the higher order ones reduce to the two- 
point one by the validity of Wick's theorem in a ballistic 
system. Alternatively, the reduced density matrix also 
completely characterizes a nonequilibrium steady state. 
The reduce density matrix, obtained by tracing over the 
bath degrees of freedom, is a better (local) quantity to 
define steady state. For the full density matrix, some sort 
of limit of bath degrees going to infinity and time going 
to infinity has to be taken in order to reach steady state, 
but such limits may not be well-defined. 

In this section, we present the method of Dhar, Saito, 
and Hanggi jl09l | who gave a procedure to compute the 
reduced density matrix in a nonequilibrium steady state. 
The starting point is an ansatz that the reduced den- 
sity matrix, although unknown, must be quadratic in the 
basic dynamic variables uc and pc of the center. For 
notational simplicity, we'll drop the subscript C in the 
following. We start by defining a vector (p = {u,p)'^, 
then we can write the reduced density matrix of the cen- 
ter as 



We note that this nonlinear self-energy is real, thus only 
shifting the frequencies of the modes, and hence it cannot 
give a finite life-time for the phonons, unable to describe 
diffusive transport. This NEGF version of "effective 
phonon" theory is closely related to the effective phonon 
theory of He et al., where the temperature-dependent 
force con stan ts are derived based on Feynman- Jensen in- 
equality [l07j . We'll call this version of mean- field theory 
as the self-consistent mean- field (SCMF). 

The current can still be calculated using Car- 
oli/Landauer formula with incorporated in C". Fig- 
ure S] shows the comparison between SCMF (solid lines) 
and master equation approac h (c rosses) in the weak 



(127) 



system-bath coupling regime [108|. Since the master 
equation approach becomes computationally very de- 
manding for the number of particles > 3, we restrict our 
comparison to one and two particle systems as shown in 
Fig. m (a) and (b), respectively. As the master equation 
formulation makes no assumptions for the strength of the 
anharmonicity, it should be considered as a numerically 
exact result, bearing in mind that the system-bath cou- 
pling is weak. Surprisingly, the SCMF approach matches 
the master equation formulation for very strong values of 
anharmonicity indicating that the self-consistency proce- 



where ^ is a matrix with twice the size of the degrees of 
freedom of the center, and the proportionality constant 
can be fixed by normalization, Tr(p) — 1. We note that 
the equilibrium statistical mechanics Gibbs distribution, 
exTp{—/3Hc), is also of this form. The nonequilibrium 
distribution introduces mixing terms between u and p. 
We determine A by matching the Green's functions. It is 
not necessary to use the complete time-displaced Green's 
functions. It is sufficient just to use the static ones, i.e., 
the Green's functions at equal time, or the covariance 
matrix 



C = (ipip^) 



{up' ) 
{pp^) 



(128) 



where the angular brackets mean trace with respect to 
the reduced density matrix p. The uu correlation can 
be computed with the greater or lesser Green's functions 
at time or using the integral of Green's function in 
frequency domain, e.g.. 



/-t-00 T 
-00 27r 



(129) 
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Using the Keldysh equation G< = G"^E<G°, S< = 
+ S^, and — —ifoJ'a, a — L,R for the two- 
lead situation, the above expression can be evaluated. 
The terms involving momenta (velocities in our con- 
vention of unit mass) can also be computed by noting 
d{u{t)u{t'Y)/dt = {p{t)u{t')) = iTidG>{t,t')/dt. Thus, 
in frequency domain, each derivative introduces an extra 
zbia; factor. 

If Lp were just ordinary vector of numbers, and the 
distribution is a general Gaussian, it is easy to see that 
A — ^C~^. This is no t tru e when if are operators. An 
important step in Ref. Il09l is the introduction of a lin- 
ear symplectic (or canonical) transform, (j)' = Sip, which 
preserves the commutation relations, and satisfies 



SJS' = J 



/ 

-/ 



(130) 



S- (C + C^) = Diag(di, • • • , dAT, di, • • • , dw), (131) 

A — S'^Diag(ai, • • • , oat, ai, • • • aN)S, (132) 

where / is the identity matrix. S is chosen such that 
the symmetrized C matrix is diagonalized with diago- 
nal elements dg repeated twice. Simultaneously, S also 
diagonalizes A with diagonal elements ag. A numerical 
procedure to do this is given in appendix of Ref. llOOi 
Since in variable (j)' the system is diagonal, the problem 
becomes equivalent to the problem of a set of decoupled 
harmonic oscillators each at a certain effective tempera- 
ture. This gives the relation 



ds = — coth(7ias 



(133) 



The exact expression for the reduced density matrix of- 
fers a good way to compare with the quantum master 
equation approach, which will be discussed in the next 
section. 



XIII. QUANTUM MASTER EQUATIONS 

In this last section before conclusion, w e tak e a look 
at the quant uinmaster equation approach [llC| to ther- 
mal transport llll - lll3 | from the point of view of NEGF. 
Although NEGF is a complete theory for answering the 
questions of thermal currents and other observables, it 
is still very difhcult to handle nonlinear systems in gen- 
eral. On the other hand, the quantum master equation 
approach handles nonlinearity with great ease: any finite 
degree center is treated the same way by expanding in 
the eigenstates of the center. But the price we have to 
pay is that we cannot treat the couplings between the 
baths and center exactly. However, if we can systemati- 
cally improve the weak-coupling approximation, then the 
master equation approach offers a great advantage. 

Ju st like NE GF, the master equations have a long his- 



Such equations can give steady state solutions for the 
dens i ty m atrix accurate to the lowest, zeroth order only 
[list Ill9l | . Some progress is made recently |l20j | in ex- 
tending the accuracy of the density matrix to second 
order by a novel analytic continuation without actually 
solving m ore complicated fourth order master equation 
[l2ll Il22| . Formally exact quantum master equation ex- 
ists either in a time-nonlocal form ^123 , 124] or time- local 
form [l25l Il26j . Here we give a transparent derivation of 
the higher order time-local master equation as well as 
energy current, using the contour order as a connection 
point to NEGF. 

We'll restrict the scope to nonequilibrium steady state 
only, although a generalization to time-dependent dy- 
namics is straightforward. There are two quantities of 
interests — one is the reduced density matrix p, and the 
other is the current II- Both of them can be treated in 
a similar way. The standard approach for obtaining the 
steady state is to evolve from the remote past product ini- 
tial state adiabatically (or even suddenly as we did in the 
earlier part of this review). In NEGF, there is no prob- 
lem with this approach since the couplings between the 
leads and the center system are handled exactly. How- 
ever, in the master equation approach (specifically, when 
the couplings themselves are treated perturbatively) , this 
adiabatic switch-on results in divergences for both the 
reduced density matrix and the current beyond the low- 
est order, generically for any initial product state pops 
where pa denotes the density matrix of the center and 
Pb — PlPr for the equilibrium baths. Unless the ini- 
tial state is carefully chosen, a steady state cannot be 
reached. 

The way to overcome this divergence is to impose a 
steady state condition for the initial state, by the re- 
quirement that the rate of change of the reduced density 
matrix p should be 0. With this, the initial state po is 
determined from a condition, rather than an initial input 
that takes any arbitrary value. Since po needs to be de- 
termined, and formally we can create a unique, invertible 
map Po to p, the problem is equivalent to determining an 
equation for p, which is commonly known as the master 
equation. 

The one-to-one map po p exists only before the adi- 
abatic limit (e — >■ 0+) is taken. Then our recipes are the 
foUowings: (1) give a formal Dyson expansion result for 
the physical quantity (O), in terms of po, where O can 
be the Hubbard operator Xmn = |'Ti)(n| (where \n) is 
the n-th eigenstate of the isolated center) or the current 
operator 1^ — pj^V^'~'uc', (2) Using (1) to obtain both p 
and dp/dt, inverting the relation from po to p and sub- 
stituting it back into the original physical quantities, we 
simultaneously obtain the equation for the current as well 
as the master equation to any desired order of accuracy 
in terms of the couplings. 

Here we introduce some notations: the total Hamilto- 
nian is H ~ h + V where h = He + Hl + is the 



tory [ll4lll5j |. The most commonly used one s keep th e decoupled system and V = u^V^'^ uc + u^V^^ uc is the 



system-bath coupling accurate in second order 



s kee p tn ( 



bath-system coupling potential. Working in the interac- 
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tion picture with respect to h, and setting the synchro- 
nization time among different pictures to 0, we have, for 
any observable 



{OH{t)) =Tr[poPB^(-oo,i)0(t)5(t,-cx))] 



Tr 



PoPbTc \ 0{t)e 



A r V(T)d7 



(134) 



where 0{t) = g^th/hQ^-ith/h -g ^j^g observable in the in- 
teraction picture, while Onit) is the same observable in 
the Heisenberg picture. The contour C now runs from 
— cx) to time of interest t. In V(t) we have included im- 
plicitly an adiabatic switch-on parameter e*^* in addition 
to the usual interaction picture form of V. The param- 
eter A = (—i/h) serves as a small expansion parameter 
in a Dyson expansion (we could also absorb a small pa- 
rameter of the coupling from V into A). We have as- 
sumed that [po,Hc] = for the validity of Eq. (|134p . 
but this is not a fundamental limitation. We can al- 
ways use p'q = Q^hto/h p^^-ihto/h ^-^j^ finite instead 
of ^0 ~^ —oo. Since po or p'q is eliminated in the end, the 
results below are independent of this assumption. We 
also note that the rate of change of O at time t is given 
by 



dt 



(0^(0) =Tr poPBTc{{dit) + 

A[O(t),y(0])e^/c^(^)'^^} 



(135) 



where [0(t),V^(t)] is the commutator of the two oper- 
ators. We'll use the symbol X to denote a matrix with 
matrix element X„„j = \n) (m|. Then the reduced density 
matrix at time t = can be computed as p = {Xj^{0)) 
and its derivative can also be similarly computed. Per- 
forming the power series expansion for the exponential, 
and noting that an odd number of bath operators or 
ur gives 0, we obtain 

p = {X^) + ^{X^V') + ^{X^V') + 0{X% (136) 

where we have introduced a short-hand notation of the 
angular brackets to mean Ti^pqPbTc ■ ■ ■] and where 
the number of contour integrals depends on the number 
n in = V{Ti)V{T2)---V{Tn). As the first term is 
explicitly po = {X^) , it is possible to invert this equation, 
to express po in terms of the final p, giving up to 6-th 
order, as 

^,{{x-v^),x-v^) + ^{{x-v'),.x-v^) 



A6 
"(2!)s 



{{{X^V^)pX^V^)X^V^)+OiX''), (137) 



where the meaning of the angular brackets are changed 
again. The brackets at the innermost level with a sub- 
script p is the same as before expect that po is re- 
placed by p. The outer slightly larger angular brack- 
ets means a trace over the density matrix of the bath, 
Pb, as well over the system with a matrix produced 
by (• • •) inside it. In addition, we still have the con- 
tour integrals. For example, let p'^^ = {X'^V'^)p, 

i.e., Pnm = TT[ppB\m){n\ J dTi J dT2TcV{Ti)V{T2)], 

then the second A"* term is ((X'^V^) „X'^V^) = 

Tr[p(^^PB\m){n\JdTiJdT2T,V{Ti)V{T2)], which is, im- 
plicitly, a linear function of p. 

The derivative of the density matrix, dp/dt, can be 
similarly expanded using Eq. ()135p . and when po is sub- 
stituted with Eq. dTF 



master equation [l21 



, we formally obtain the time-local 
122|, up to 6-order, as 



dp 



-^^■p + A^ ([X^, V]V), + -([X^, V]V^)p 
-^{{X^V^)^[X\V]V) + ^{[X^^V]V'), 

-^{ix^v')^[x^.y]V) + 

A6 
2! 2! 

2! 3! 



+ ^,{{^X^^^)pX^V^)[X^., V]V) + 

-{{X^V^)p[X^, V]V^) + 0{\^) = 0. (138) 



The meaning of the first term is — z 



E„-E„ 



if we 



write out explicitly in matrix element form, where £"„ is 
the eigen-energy of the n-th state of the isolated, arbi- 
trary nonlinear center of He- The meaning of the two 
types of angular brackets remains the same. The time 
arguments for are at t = 0, while all the other 

Vs have dummy contour time argument and need to 
be integrated out. The p dependence is in the angular 
brackets (• • •)p = Tr [pps^c J^dr ■ ■ ■]. After performing 
the trace and contour integrals, we obtain explicitly the 
equation for p. If we truncate the series to second or- 
der in A, we g et the standard Redfield quantum master 
equation |ll6{ . 

The current can be treated in a similar fashion. In 
fact, the mathematical structure of the current is the 
same as that of master equation, except that we just need 
to replace the commutator, [X'^,V], hy V = pJ^V^'^uc 
where the left-sided dot on V indicates that the time 
derivative is performed only for the left lead, then we 
can write, for the quantity / = XI l, as 



(/h> = A(V> =A(Ve^/c^("^''") 



(139) 



X4 

3! ^ 
A4 



5! 



= Vy>p + ^-^{VV')p ^{{X''v')pVV) + 0(Ae). 



The last line is due to Eq. (jl37l) where po is written in 
terms of p. A divergence appearing in the second term 
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<VV> = -Q-iC^ + <-522i-o- 

k I j k 

(i) (ii) 

(1) (2) (ST^ 



+ -° 




^«xV'>vv> 



(10) (11) (12) 

(a) (b) (c) 



(d) (e) (f) 

FIG. 5. Diagrams representing the terms in the current. The 
graphs (1), (4), (5), (7), (8), (10) and (a)-(f) have divergent 
terms of the form cx 1/e. The Feynman rules are discussed in 
the text. 



gets cancelled explicitly by the third term. This approach 
solved the divergence problem which has been puzzling 
us for a while. The above result is a general izati on of the 
seco nd o rder result (the first term) in Ref. Il27l see also 
Ref. HH. 

It is possible to represent various terms for the expres- 
sion of current in terms of diagrams after unraveling the 
contour time into normal time with time order or anti- 
time order. The diagrams are shown for the various terms 
in Fig. [5] The Feynman rules for the diagrams are as fol- 
lows: 1) each dot is associated with a time tj and the 
system operator; each line segment between the dots has 
a system state label; the matrix element of the operator 
is (fc|wclOe*'^''''*'■^'*^ where we define A^; = {Ek-Ei)/h. 
2) The dots are connected by the phonon lines, represent- 
ing the contour function C(ti, T2) — ihT,{Ti, T2), tj = tj, 
in all possible ways. 3) The open dot has a fixed time of 0, 
and associated C has a time derivative. 4) The reduced 
density matrix p is represented as a square box. 5) The 
direction of the arrow represents ordering; right pointing 
arrows are for anti-time order (lower branch), and left 
pointing arrows time order (upper branch). The hori- 
zontal line represents the trace over the system states. 
With these rules, e.g., the first two diagrams, (i), (ii), 
and diagram (3) are (assuming center has only one de- 
gree of freedom and Ski — {k\uc\l}, for multiple degrees 
of freedom, C becomes a matrix and S received a vector 
index) 



(ii) = ^ / dti pkiSijSjke'^^''*'+'*' 0(0, t+), (141) 
kij •'-'^ 

(3) = ^ PkiSipSpgSqjSjk I dti I dt2 I dt3 (142) 

klpqj 



— 00 J —oc J —CO 



We intend to give a full account of this method with 
numerical application elsewhere [129]. 



XIV. CONCLUSION 

In this review, we gave simple examples of Green's 
functions for harmonic oscillator as a starting point. We 
have tried to emphasize the contour ordered functions 
as the basic language of NEGF. The functions are de- 
fined on a finite segment from the initial time to the 
current time of interest. This formulation makes the the- 
ory work equally well for steady state and transient time 
development. A number of applications explore this fea- 
ture, noticeably the transient problem of full counting 
statistics. The standard topics of NEGF, the equation of 
motion method and Feynman-diagrammatic expansion, 
were briefly discussed. Some recent developments are re- 
viewed, such as the new formulas for coupled left-right 
leads. In the treatment of nonlinear systems, we draw at- 
tention to the self-consistent mean field theory, which was 
shown to give surprisingly accurate result for the current 
for small systems. We attempted to make connections 
between NEGF and master equation. The master equa- 
tion approach focuses on the reduced density matrix. An 
exact expression for the density matrix for ballistic sys- 
tem is reviewed. The last section is a bit off the main 
line of this review. There, we gave a very transparent 
derivation of the higher order (time-local) quantum mas- 
ter equation and a suggestion on how higher order current 
can be computed. This last part is new, to our knowl- 
edge. 
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E 



dh PkiSi,Sjke''^'^''+'''C{0,t^ 



(140) 
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